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3-D  Numerical  Simulations  of  Biofilm  Dynamics  with  Quorum 

Sensing  in  a  Flow  Cell 

Jia  Zhao*and  Qi  Wang* 


Abstract 

Biofilms  are  microorganisms,  where  bacteria  are  embedded  in  networks  formed  by  bacteria-produced 
expopolysacharrides,  also  known  as  extracellular  polymeric  substances  (EPS).  Bacteria  in  biofilms  com¬ 
municate  and  cooperate  with  each  other  by  sensing  the  density  of  signalling  molecules  that  bacteria 
excrete.  This  phenomenon  is  known  as  quorum  sensing  (QS).  In  this  paper,  we  develop  a  hydrody¬ 
namic  model  for  biofilms  to  include  quorum  sensing  in  biofilm  dynamics,  extending  previous  models 
for  biofilms.  In  the  new  model,  we  classify  bacteria  into  down-regulated  and  up-regulated  quorum  sens¬ 
ing  cells,  as  well  as  non  quorum  sensing  cells  based  on  their  quorum  sensing  ability.  We  then  develop 
a  numerical  simulation  tool  for  studying  biofilms  under  quorum  sensing  in  an  aqueous  environment  by 
solving  the  hydrodynamic  model  in  full  3D  in  space  and  time.  Our  numerical  results  show  that  this  model 
captures  key  features  of  quorum  sensing  regulation  that  are  responsible  for  the  development  of  heteroge¬ 
neous  biofilm  structures.  For  instance,  numerical  simulations  with  the  model  show  that  quorum  sensing 
is  beneficial  for  the  biofilm  development  in  a  long  run  by  building  a  robust  EPS  population  to  protect 
the  biofilm,  but  maybe  of  little  benefit  in  some  short  time  frame  because  it  slows  down  the  development 
of  biofilm  colonies.  In  addition,  numerical  simulations  demonstrate  that  quorum  sensing  induction  in 
biofilms  is  sensitive  to  the  hydrodynamic  stress  and  competes  directly  with  the  nutrient  supply.  When 
nutrient  supply  is  strong  in  the  entire  domain,  QS  induction  produce  more  QS  up-regulated  bacteria  to 
grow  robust  biofilm  colonies  upstream.  On  the  other  hand,  when  nutrient  supply  is  weak,  QS  induction 
is  more  prominent  downstream.  Hydrodynamic  stress  also  alters  morphology  of  the  biofilm  when  the 
flow  is  strong  in  the  domain. 


1  Introduction 

Bacteria  are  ubiquitous  in  nature  as  well  as  in  our  daily  life.  In  general,  bacteria  usually  do  not  exist 
as  isolated  cells  but  rather  live  in  organized  communities.  These  bacterial  communities,  where  bacteria 
are  less  motile  and  are  glued  together  by  exopopolysacharrides  or  extra  cellular  substances  (EPS),  are  called 
biofilms.  It’s  commonly  perceived  by  the  medical  community  that  biofilms  are  responsible  for  many  diseases 
or  ailments  associated  with  chronic  infections,  supported  by  the  survey  data  that  biofilms  are  present  on 
the  removed  tissue  of  80%  of  patients  undergoing  surgery  for  chronic  sinusitis  [?].  Unlike  a  planktonic 
(free- swimming)  bacterium,  bacteria  in  biofilms  appear  to  be  more  adaptive  to  the  stressful  environment,  for 
instance,  they  are  more  tolerable  to  antimicrobial  agents  and  hydrodynamic  stress.  As  the  result,  biofilms  are 
always  hard  to  be  eradicated  by  standard  antimicrobial  treatment  [?],  which  perhaps  explains  the  tendency 
for  relapse  of  chronic  diseases  or  ailments  caused  by  biofilms. 
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One  of  the  main  features  of  biofilms,  which  make  bacteria  in  biofilms  functioning  quite  different  from 
planktonic  bacteria,  is  that  bacteria  living  in  biofilms  are  observed  to  communicate,  and,  further  more,  co¬ 
operate  with  each  other  by  secreting  solution-dissolvable  signalling  molecules  known  as  autoinducers  [?]. 
For  gram-negative  bacteria,  this  particular  molecule  is  acyl  homoserine  lactones  (AHLs)  [?].  Experimental 
results  have  shown  that  bacteria  can  sense  the  concentration  of  the  signaling  molecules  secreted  by  bac¬ 
teria  in  the  surrounding  and  function  cooperatively  to  accomplish  certain  tasks  when  the  concentration  of 
antoinducers  reaches  a  threshold  value  [?].  This  phenomenon  is  commonly  known  as  quorum  sensing  (QS). 
The  signalling  molecule  or  the  autoinducer  is  therefore  called  the  quorum  sensing  molecule.  This  cell-cell 
communication  system  (via  quorum  sensing)  has  been  observed  in  many  bacterial  systems,  such  as  vibrio 
fischeri  [?],  where  quorum  sensing  was  first  observed,  Pseudomonas  aeruginosa,  Escherichia  coli  and  so  on. 
This  approach  allows  unicellular  microorganism  to  perform  behavior  like  multicellular  organism. 

Quorum  sensing  has  been  shown  to  be  responsible  for  mediating  a  variety  of  social  activities  in  biofilms, 
which  include  the  secretion  of  diverse  byproducts,  biofilm  growth  [?],  swarming  motility  and  virulence  gene 
expressions  [?].  For  an  overview  of  quorum  sensing  mechanism  in  biofilm  formation  and  function  and  up- 
to-date  advances,  readers  are  referred  to  the  work  of  Miller’s  [?]  and  Waters’s  [?].  In  addition,  some  non¬ 
signaling  properties  of  quorum  sensing  molecules  have  also  been  reported  in  the  literature  [?].  Especially, 
the  phenomenon  of  quorum  sensing  regulating  expopolysacharride  production  during  biofilm  formation  has 
been  reported  widely  in  the  literature  [?,  ?,  ?]. 

EPS  can  be  treated  as  the  so-called  ‘public  good’  since  its  a  shared  resource  that  each  individual  bac¬ 
terium  can  take  advantage  of,  even  if  it  is  produced  by  other  bacteria.  However,  how  quorum  sensing 
regulation  affect  the  EPS  production  is  not  clearly  resolved.  Besides,  in  some  bacterial  systems,  like  Vibrio 
cholerae,  the  EPS  production  is  suppressed  when  the  density  of  autoinducer  reaches  its  threshold  value  [?]; 
whereas  in  Pseudomonas  aeruginosa,  the  EPS  production  would  be  facilitated  once  the  quorum  sensing  is 
activated  [?].  In  addition,  as  the  phenomenon  involves  transport  of  signaling  molecules,  the  heterogenous 
structure  in  biofilm  colonies  and  the  solvent-biomass  interaction  can  not  be  ignored,  as  they  can  affect  the 
convection  and  diffusion  rates  of  the  signaling  molecules  transported  as  well  as  the  morphology  of  biofilm 
colonies.  This  motivates  us  to  develop  a  full  3D  mathematical  model  to  study  how  quorum  sensing  regulates 
biofilm  formation  and  development  as  well  as  the  pros  and  cons  of  quorum  sensing  in  an  aqueous  environ¬ 
ment  where  hydrodynamic  interaction  between  the  various  biofilm  components  and  the  ambient  fluid  flow 
is  important. 

The  induction  of  quorum  sensing  is  usually  observed  at  high  bacterial  cell  volume  fractions,  where 
the  concentration  of  QS  molecules  is  accumulated  to  a  certain  threshold  value.  However,  a  high  bacterial 
cell  density  does  not  necessarily  induce  quorum  sensing  since  many  other  factors  have  been  observed  to 
contribute  to  this  induction  such  as  nutrient  supplies,  hydrodynamic  shear  stress  [?]  and  the  PH  environment. 
Intuitively,  a  hydrodynamic  flow  may  have  a  diluting  effect  on  the  signalling  molecules,  which  then  lead  to 
a  delay  or  even  failure  of  QS  induction  by  preventing  QS  molecules,  like  AHLs,  from  reaching  a  certain 
threshold  value.  For  instance,  in  [?],  the  author  claims  that  the  amount  of  biomass  required  for  the  full  QS 
induction  of  the  population  increases  as  the  flow  rate  increases.  On  the  other  hand,  coming  with  the  inflow 
solvent,  more  essential  substances,  such  as  nutrient  which  fuels  the  production  of  QS  molecules  by  bacteria, 
would  be  fresh  in;  thus  hydrodynamic  flow  may  lead  to  higher  production  of  the  signalling  molecules  and 
thereby  facilitate  QS  induction.  Little  work  has  been  done  to  analyze  the  correlations,  not  to  mention  taking 
the  heterogenous  structure  factor  of  biofilms  into  account.  This  adequately  serves  as  another  motivation  for 
us  to  develop  a  full  3D  hydrodynamic  biofilm  model  to  study  the  hydrodynamic  effect  on  quorum  sensing. 

In  the  literature,  related  with  quorum  sensing,  some  other  factors  have  also  been  singled  out  for  possible 
impact  to  biofilm  development.  For  instance,  recently,  it  has  been  verified  that,  except  for  cooperation, 
there  probably  exist  cheaters  in  the  bacterial  population  which  exploit  the  signaling  molecules  produced 
by  others  [?].  Vermant’s  group  observed  that  quorum  sensing  coordinated  secretion  of  rhamnolipid  acts  as 
a  surfactant,  which  leads  to  finger  pattern  formation  in  swimming  bacteria  aggregates  [?].  With  the  ever- 
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increasing  knowledge  about  quorum  sensing,  researchers  are  working  on  developing  potential  applications 
by  taking  advantage  of  this  cell-cell  communication  mechanism.  For  an  overview  of  applications  of  quorum 
sensing  in  biotechnology,  readers  are  referred  to  [?].  One  promising  application  is  that  it  is  probably  more 
efficient  to  use  anti-quorum  sensing  treatment  than  antibiotic  agents  in  treating  some  biofilm  infections, 
which  is  likely  to  pose  a  selective  pressure  for  the  development  of  resistant  mutants  [?]. 

Inspired  by  experimental  findings,  researchers  have  come  up  with  some  mathematical  models  to  study 
biofilm  formation  and  function,  as  well  as  related  issues,  such  as  the  mechanism  of  quorum  sensing  and 
its  contributions  to  biofilm  structures  and  activities.  Two  review  papers  [?,  ?]  have  given  an  overview  of 
current  advances  in  mathematical  modeling  of  biofilms.  Concerning  coupling  biofilm  growth  with  quorum 
sensing  features,  several  mathematical  models  have  been  developed,  some  of  which  are  cell-based  [?],  some 
are  spatially  homogeneous  consisting  of  ordinary  different  equations,  which  are  derived  based  on  the  mass- 
action  principle  [?]  while  others  are  models  consisting  of  partial  differential  equations.  In  particular,  John 
Ward  derived  a  series  of  models  for  analyzing  quorum  sensing  mechanisms  in  one  space  dimension  [?,  ?]. 
Hermann  derived  a  2D  partial  differential  equation  model  coupled  with  the  Stokes  flow  [?]  and  later  he 
extended  this  model  [?]  and  proved  its  wellposedness  [?].  In  [?],  the  author  proposed  a  2D  quorum  sensing 
model  and  studied  the  shape  of  biofilms  at  the  onset  of  the  hydrodynamic  shear.  However,  there  is  very  little 
work  on  3D  hydrodynamics  of  biofilms  and  the  ambient  fluid  environment  available  in  the  literature,  which 
takes  into  account  the  spatial-temporal  heterogeneous  structure  of  biofilms  as  well  as  the  hydrodynamic  flow 
effects.  This  is  perhaps  mainly  because  of  the  mathematical  complexity  and  computational  challenges  in  re¬ 
solving  biofilms  which  are  a  truly  3D  microorganism  with  highly  heterogeneous  spatial-temporal  structures 
and  complex  dynamics. 

In  this  paper,  we  develop  a  new  full  3D  hydrodynamic  model  for  the  mixture  system  consisting  of 
biofilms  and  the  ambient  fluid,  extending  our  previous  phase-field  model  for  biofilms  of  two  and  three  ef¬ 
fective  components  [?,  ?].  This  model  couples  hydrodynamics  of  the  biofilm-solvent  mixture  to  the  quorum 
sensing  activity,  as  well  as  dynamics  of  nutrient  and  antimicrobial  agents.  Using  this  model,  we  aim  to  in¬ 
vestigate  biofilm  formation  regulated  by  quorum  sensing  in  an  aqueous  environment  under  hydrodynamical 
flows  and  possibly  antimicrobial  treatment.  A  second  order  numerical  scheme  in  both  time  and  space  based 
on  finite  difference  methods  is  devised  to  solve  the  governing  system  of  partial  differential  equations  in  the 
model.  The  numerical  solver  for  the  discretized  equations  is  implemented  on  graphical  processing  units 
(GPUs)  in  a  3D  cubic  geometry  in  space  and  time. 

The  3D  numerical  simulation  tool  is  put  in  use  to  investigate  the  role  of  quorum  sensing  to  the  dynamics 
of  various  biomass  components  as  well  as  to  the  entire  biofilm  colony.  Given  the  complexity  of  the  problem, 
we  select  a  couple  of  important  case  studies  to  present  in  this  paper.  Firstly,  we  benchmark  the  model  on  two 
features  of  quorum  sensing:  the  sensitivity  of  biofilm  development  to  the  threshold  value  of  quorum  sensing 
induction  and  the  effect  of  the  growth  rate  of  the  quorum  sensing  down  regulated  bacteria.  The  second  case 
is  chosen  because  after  having  conducted  many  numerical  experiments  on  model  parameters,  we  notice  that 
the  model  is  very  sensitive  to  the  growth  rate  of  the  QS  down-regulated  bacteria.  Then,  we  apply  the  full  3D 
simulation  tool  to  investigate  the  interaction  between  the  hydrodynamics  and  quorum  sensing  regulation. 
Finally,  we  examine  the  impact  of  hydrodynamic  stresses  on  biofilm  morphology. 

The  rest  of  this  paper  is  organized  into  four  sections.  In  the  second  section,  we  present  the  detail  of  the 
derivation  of  the  hydrodynamic  phase  field  model  for  biofilms.  Then,  we  design  an  efficient  numerical  solver 
for  the  governing  partial  differential  equations  in  the  model  using  semi-implicit  finite  difference  schemes 
in  the  third  section.  In  the  fourth  section,  the  numerical  results  and  discussions  are  presented.  In  the  last 
section,  we  conclude  the  study. 
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2  Mathematical  model 


In  this  section,  we  extend  a  previously  developed  hydrodynamic  model  for  biofilms  [?]  to  include  the 
effect  of  quorum  sensing  among  the  various  bacterial  phenotypes  classified  based  on  their  quorum  sensing 
abilities.  The  biofilm,  consisting  of  bacteria,  EPS  and  solvent,  is  coarse-grained  into  a  single  fluid  model 
with  multi-components,  in  which  different  types  of  bacteria,  EPS,  as  well  as  solvent  are  modeled  effec¬ 
tively  as  separate  fluid  components.  Quorum  sensing  molecules,  nutrient  and  antibiotic  agents  are  treated  as 
phantom  materials  which  exhibit  chemical  reactive  effects  without  their  mass  and  volume  being  considered. 
Treating  a  small  amount  of  substances  of  such  small  molecules  as  phantom  materials  is  an  effective  approx¬ 
imation  to  the  complex  mixture  system  given  the  minute  contribution  from  these  materials  to  the  total  mass 
and  volume  of  the  material  system. 

2.1  Basic  notations 

We  study  biofilms  in  a  cubic  domain  fi,  where  Lx,Ly,  Lz  represent  the  length  in  the  x,y,z  direction, 
respectively,  i.e., 

fi=[0  Lx\  x  [0  Ly\  x  [0  Lz\. 

Following  our  previous  approach  [?,  ?],  we  denote  the  volume  fraction  of  the  effective  solvent  and  biomass 
as  (j)s  and  0n,  respectively.  In  biomass,  we  further  divide  it  into  EPS  and  bacteria,  whose  volume  fractions 
are  denoted  by  <fip  and  05,  respectively.  By  definition,  it  follows  that  0n  =  0P  +  05.  The  incompressibility 
condition  of  the  mixture  system  requires 


05  +  4>p  +  06  —  1*  (2.1) 

Let  pb,  v&,  vp  and  ps ,  v5  be  the  density  and  velocity  of  the  bacteria,  EPS  and  solvent,  respectively.  Then 
the  volume-averaged  density  and  velocity  for  the  mixture  is  denoted  as 


P  —  0 bpb  T"  (frpPp  T"  0 5 Ps  5  V  —  06  V6  “1“  0£>VJ9  05^5*  (2.2) 

In  the  biofilm  model,  we  classify  the  bacteria  into  four  different  phenotypes  according  to  their  responses 
to  quorum  sensing  (QS)  molecules  [?,  ?]:  QS  down-regulated  bacteria  ( QS~ ),  QS  up-regulated  bacteria 
(QS"+),  non-QS  bacteria  (non-QS"),  as  well  as  the  dead  bacteria,  whose  volume  fractions  are  denoted  by 
051 ,  052,  063  and  064,  respectively.  Hence,  the  volume  fraction  of  the  bacteria  is  given  by  the  sum  of  the 
volume  fractions  of  the  four  types: 


06  —  061  +  062  +  063  +  064-  (2.3) 

In  biofilms,  QS~  are  the  bacteria,  which  can  convert  into  QS+  once  the  concentration  of  QS  molecules 
reaches  a  threshold  value;  whereas  QS+  are  the  bacteria,  which  can  convert  into  QS~  when  the  concen¬ 
tration  of  QS  molecules  is  lower  than  the  threshold.  In  another  word,  QS~  bacteria  belong  to  the  wild 
phenotype  of  the  bacteria  and  QS+  bacteria  represent  the  phenotype  who  has  bonded  QS  molecules  and 
certain  activated  gene  transcription  for  enhanced  EPS  production  and  therefore  reacts  differently  from  its 
original  phenotype.  Non-QS"  bacteria  represent  those  bacteria  who  are  not  affected  by  the  concentration  of 
quorum  sensing  molecules  but  may  be  regulated  by  other  proteins  in  the  biofilm.  These  non-quorum  sensing 
bacteria  could  be  regarded  a  phenotype,  which  do  not  prefer  a  social  relation  with  other  bacteria,  or  they 
may  also  be  treated  as  a  different  type  of  bacteria,  such  as  invaders.  Dead  bacteria  are  those  who  lose  their 
vitality  but  are  still  attached  in  the  biofilm. 

In  addition  to  the  above  bacterial  phenotypes,  we  also  keep  track  of  the  functional  components  whose 
mass  and  volume  are  omitted  in  the  model,  but  their  chemical  effects  are  retained  as  their  mass  fractions  are 
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Non-QS  Bacteria 


Dead  Bacateria 
EPS 


0  Antibiotics 
^  QS  Molecue 
if  Nutrient 


Figure  1:  A  Schematic  cartoon  for  a  biofilm  colony  with  multiple  components. 


much  smaller  than  those  of  the  biomass  and  solvent.  These  are  nutrients,  antimicrobial  agents  and  quorum 
sensing  molecules,  whose  concentration  are  denoted  by  c,  d  and  A,  respectively.  To  help  the  reader  better 
grasp  the  components  in  our  model,  we  show  a  schematic  cartoon  for  the  biofilm  colony  in  Figure  ??. 

We  make  a  simplifying  assumption  here  that  all  bacterial  cells  and  EPS  that  they  produce  interact  with 
the  solvent  equally  in  terms  of  thermodynamics  (entropic  motion  and  mixing).  Thus  the  effective  free  energy 
density  functional  of  the  biological  system  can  be  proposed  as  follows  using  the  extended  Flory-Huggin’s 
mixing  free  energy  density  with  an  extra  conformational  entropy  term  [?] 


f  =  kBT(11  ||V</>n||2  +  72 


log  (j)n  +  (ps  log  (ps  +  X(Pn(ps 


(2.4) 


where  is  the  boltzmann  constant,  T  is  the  absolute  temperature,  N  is  an  effective  polymerization  index 
for  the  non-solvent  biomass  components,  71  parametrizes  the  strength  of  the  conformational  entropy,  72 
parameterizes  the  strength  of  the  bulk  interaction  potential,  x  is  the  mixing  coefficient.  A  more  refined 
potential  could  be 


f  =  kBT  7l||V<)>n||2  +  72 


loS  (f>b+  lo§  fip  +  (t>s  loS  (t>s  +  X(pn(ps 


Nb 


(2.5) 


where  N &  and  Np  are  two  effective  polymerization  index  for  the  bacteria  and  EPS,  respectively.  In  the 
problems  we  are  interested  in,  the  difference  between  the  two  potential  is  negligible.  So,  we  adopt  the 
simplified  one  in  (??). 


2.2  Transport  equations  for  biomass  components 

With  the  proposed  thermodynamic  free  energy,  the  transport  equation  for  each  bacterial  cell  type  is 
assumed  to  be  governed  by  a  reactive  Cahn-Hilliard  type  equation: 

+  V  •  (v(pbj)  =  V  ■  (XjfpbjV fij)  +  gb],  j  =  1,2, 3, 4,  (2.6) 
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where  X3  is  the  mobility  parameter  for  bacterial  type  j,  fi3  —  is  the  chemical  potential  with  respect  to 
bacterial  type  j ,  and  g\)3  is  the  reactive  rate  for  cell  type  j. 

The  population  dynamics  of  these  cell  types  include  regulated  growth  due  to  nutrient  and  decay  due 
to  antimicrobial  agents  as  well  as  natural  causes.  Between  the  quorum  sensing  cells,  there  exists  the  up 
and  down  regulation  coordinated  by  the  concentration  of  the  quorum  sensing  molecule  AHL.  The  reactive 
kinetics  for  all  cell  types  are  given  by 

9bi  =  (■ Cbi<f>bi  +  t3Cb2(f>b2)k^rchi^  ~  ^77)  -  T<t>b2  ~  ri4>bi  ~ 

9b2  =  (1  -  -  £-)  +  ftgLfa  - 

Sb 3  =  ~  “  r3^&3  “ 

9b4  =  Si=l  ri(t>bi  +  Sj=l  k^+d&bi  ~  r4^4i 

(2.7) 

where  the  C bi  are  the  growth  rates,  (3  represents  the  percentage  of  QS+  bacteria  that  produces  QS~  bacteria, 
Cdi  is  the  decay  rate  due  to  the  antimicrobial  treatment,  r{  is  the  natural  death  rate,  k kdi  are  the  half 
saturation  constants  in  the  respective  Monod  model  for  each  bacterial  phenotype,  respectively,  where  i  = 
1,  2,  3.  Here  r 4  is  the  natural  conversion  rate  from  dead  bacteria  into  solvent.  In  addition,  ru  Trf^^n  is  the 
conversion  rate  from  QS~  to  QS+ ,  and  r<i\  rJ+An  conversi°n  rate  from  QS+  to  QS~ ,  <j)b,max  is  the 
carrying  capacity  for  the  bacterial  volume  fraction,  r  is  the  threshold  AHL  concentration  locally  required  for 
quorum  sensing  induction  to  occur.  The  exponent  n  describes  the  degree  of  polymerization  in  the  synthesis 
of  AHL,  where  we  choose  n  =  2  in  the  study.  When  A  >  r,  quorum  sensing  induction  takes  place  to 
increase  the  up-regulated  cell  population  and  decrease  the  down-regulated  cell  population.  The  non-QS 
cells  are  not  affected  by  the  quorum  sensing  molecules  by  definition. 

The  transport  equation  for  EPS  which  is  the  product  of  the  live  bacteria,  is  proposed  as  a  reactive 
Cahn-Hilliard  type  equation  as  well  and  is  given  by, 

-7^  b  V  •  (v0p)  =  V  •  (A p(j)pS7 ftp)  +  gp,  (2.8) 

where  Xp  is  the  mobility  parameter,  gp  —  is  the  chemical  potential  with  respect  to  (j)p ,  and  gp  is  the 
reactive  term  proposed  as  follows 


9p  — 
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CpiC 
kpi  T  c 


i= 1 


1  - 


up,  max 


+  ar4(/>64  -  rp^cP 


kcp  +  c 


CdPd 
k dP  T-  d 


(2.9) 


Here  CP3  is  the  EPS  production  rate  with  respect  to  the  jth  type  cell,  j  =  1,2,3,  facilitated  by  the  nutrient. 
0p,max  gives  the  upper-bound  for  the  EPS  production.  We  assume  a  part  of  dead  bacterial  population  is 
converted  into  EPS  while  others  into  the  solvent,  where  a  is  the  efficiency  for  the  conversion,  ranging 
between  0  and  1.  We  also  assume  that  some  EPS  can  be  dissolved  into  the  solvent  naturally  and  accelerated 
by  the  drug  with  the  rate  given  by  rp  and  Cdp,  respectively.  kpi,  kap ,  kcp  and  kdp  are  the  half  saturation 
constants  for  the  Monod  models  adopted. 

Because  of  the  simplification  assumption  we  made  on  the  bacteria  and  EPS  in  the  free  energy,  the 
chemical  potential  for  each  species  is  identical  in  this  model.  So  the  interfacial  force  due  to  the  chemical 
potential  to  the  momentum  transport  is  given  by 

4 

07T,V/i77,  =  ^  ^  (j)bj  V /%  •  (2.10) 

3= 1 
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With  the  proposed  transport  equations,  we  identify  the  transport  velocity  for  bacteria,  EPS  and  solvent, 
respectively,  as  follows: 


1 


where  /in  =  ^  is  the  chemical  potential  for  (j)n.  Notice  that  the  transport  equation  for  (j)n  can  be  derived 
by  summing  up  the  transport  equation  for  each  biomass  component, 


[v  -  (f>pVp  -  <f>bVb].  (2.11) 


Vp  —  v  APV  fir 


2=1 


9(f)n 

dt 
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+  V  •  (v0n)  =  V  •  (A 0nV/in)  + 

2  =  1 


(2.12) 


where  A</>n  =  ^2^2  +  Ap</>p.  Since  we  need  to  keep  track  of  each  individual  species  in  the  model,  this 

equation  is  not  used  in  this  model. 


2.3  Transport  equations  for  the  functional  components 

There  are  three  functional  components  in  our  biofilm  model,  namely,  the  nutrient,  the  antimicrobial 
agent,  and  the  quorum  sensing  molecule,  whose  concentrations  are  denoted  by  c,  d  and  A,  respectively. 
These  molecules  all  have  small  molecular  weight  compared  with  the  components  of  the  biomass,  thus  their 
molecular  mass  and  volume  are  neglected  in  this  model  for  simplicity.  Instead  of  tracking  their  volumes, 
their  concentrations  are  traced,  which  are  governed  by  traditional  convection-diffusion-reaction  equations. 
These  molecules  are  dissolved  in  solvent  so  that  they  are  transported  with  the  solvent. 

Specifically,  the  transport  equation  for  the  concentration  of  the  quorum  sensing  molecule  AHL  (A)  is 
given  by 

+  V  •  (vs0sA)  =  V  •  (<PsDaVA)  +  ga,  (2.13) 

where  Da  is  the  diffusion  coefficient  and  ga  is  the  reactive  term,  proposed  as 


9a  =  ~raA  + 


1  - 


A 


Cgd 
kda  +  d 


A. 


(2.14) 


Here  ra  is  the  decay  rate  of  the  AHL  concentration  due  to  binding  with  QS  bacterial  cells,  ra\  and  ra 2  are  the 
growth  rate  due  to  the  down-regulated  and  up-regulated  quorum  sensing  cells,  respectively.  (  \  and  Ca 2  are 
growth  rates,  which  take  into  account  the  effects  that  the  concentration  of  AHL  would  affect  the  productivity 
of  bacteria  for  AHL  positively.  Here  Amax  gives  the  upper-bound  for  the  concentration  of  the  auto-inducers. 
Ca  is  the  decay  rate  of  AHL  due  to  antimicrobial  agents  and  //,/,,  is  the  half  saturation  constant.  The  diffusion 
coefficient  is  proposed  by  following  the  Hinson  model  [?], 


Da  =  Dao 


1  4*b  $s 

1  +  \4>b  <j)s  +  ^  ’ 


(2.15) 


where  the  first  term  on  the  right  hand  side  represents  the  reduction  effect  due  to  the  presence  of  bacteria 
while  the  second  term  is  an  empirical  fitting  accounting  for  the  reduction  effect  due  to  the  presence  of  EPS. 
Here  Dpr  is  a  model  parameter  fitted  experimentally. 

The  transport  equation  for  the  nutrient  is  proposed  as  follows 


d(j)sc 

dt 


+  V  •  (0svsc) 


V  •  (Ds<j)s\7c)  +  gc , 


(2.16) 
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where  Ds  is  the  diffusion  coefficient  for  nutrient  transport  and  the  reactive  rate  is  given  by 


9c  = 


3 

rbic 

“  kbi  +  C1 
2=1 


E 

2=1 


kpi  +  c 


(2.17) 


The  first  term  in  gc  represents  the  nutrient  consumption  due  to  bacteria  metabolism  and  the  second  term 
represents  the  nutrient  consumption  due  to  the  EPS  production. 

For  the  diffusion  rate  of  the  nutrient  (Dc),  which  is  primarily  oxygen  in  this  model,  Stewart  [?]  suggests 
the  diffusion  rate  in  biofilms  is  around  75%  than  that  in  the  pure  solvent.  Since  the  molecular  weight  of 
oxygen  is  relatively  small,  it  is  supposed  to  penetrate  EPS  and  the  membrane  of  cells  equally.  So  we  set  the 
diffusion  rate  of  oxygen  in  biofilms  as  follows 


Dc 


1  (06  0p) 

1  +  ^(06  +  0p) 


(2.18) 


where  the  second  term  in  the  expression  represents  the  reduction  in  diffusion  due  to  the  presence  of  the 
biomass. 

Finally,  the  transport  equation  for  antibiotic  agents  is  proposed  as  follows 


+  v  •  (<j>svsd)  =  V  •  {Dd<t>sVd)  +  gd, 

where  Dd  is  the  diffusion  rate  of  antimicrobial  agents  and 


(2.19) 


9d 


d 

kd  +  d 
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Cy  ^  ^di&bi  H-  T'dpfip)  T’dd', 
2=1 


(2.20) 


is  the  reactive  term.  The  antibiotic  is  assumed  to  be  consumed  by  all  bacterial 
rates.  Here  represents  the  natural  decay  rate  of  antimicrobial  agents.  For  the 
adopt  the  Hinson  model  [?] 


Dd  = 


Ddo 


1  06  0s 

1  +  (t>s  +  ^ 


cells  and  EPS  at  different 
diffusion  rate  Dd ,  we  also 

(2.21) 


In  the  derivation  of  these  transport  equations,  we  try  to  include  as  many  biological  and  physical  factors 
as  possible  which  results  in  a  large  set  of  model  parameters.  Some  of  them  can  be  measured  or  searched 
from  the  existing  literature,  while  others  are  obtained  either  using  our  best  guesses  or  calibrated  against 
some  experiments.  Nevertheless,  these  effects  are  quantitatively  present  in  the  biofilm  system. 


2.4  Momentum  and  Continuity  Equation 

In  order  to  couple  the  biofilm  components  to  hydrodynamics,  we  need  the  governing  equation  for  the 
averaged  velocity  v.  As  we  mentioned  above,  we  treat  this  biological  system  as  one  fluid  with  multiple 
components.  By  assuming  the  fluid  mixture  as  solenoidal,  we  impose  the  balance  of  linear  momentum  and 
conservation  of  mass  as  follows 

=  —Vp  +  V  •  T  -  0nV/in  (2.22) 

=  0,  (2.23) 

where  p  is  the  averaged  density,  p  is  the  hydrostatic  pressure,  r  is  the  extra  stress  and  the  last  term  is  the 
interfacial  force  due  to  the  inhomogeneity  of  biomass  distribution,  derived  from  the  virtual  work  principle. 


P  l  —  +  v  •  Vv 


dt 


Vv 
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In  this  paper,  we  propose  the  mixture  as  an  extended  Newtonian  fluid,  then  the  extra  stress  tensor  is  given 
by 

r  =  2rjb(pb'Db  +  +  2r/s(/>sDs,  (2.24) 

where  775,  rjp  and  ps  are  the  viscosity  of  bacteria,  EPS  and  solvent,  and  D5,  Dp  and  Ds  are  the  rate  of 
deformation  tensor,  correspondingly, 

Df>  =  ^(Vv&  +  VvD>  dp  =  ^(Vvp  +  VvJ)>  =  ^(Vvs  +  VvJ),  (2.25) 

where  v&,  vp,  v5  are  the  velocities  for  bacteria,  EPS  and  solvent.  Their  expressions  are  given  in  equation 

(??). 


2.5  Boundary  Conditions 

In  order  to  simulate  biofilms  in  a  fixed  domain,  we  need  to  impose  boundary  conditions  for  the  governing 
partial  differential  equations,  which  mimic  common  situations  for  biofilms  being  cultured  and  observed  in 
vitro.  In  this  paper,  we  are  interested  in  simulating  biofilm  flows  in  two  special  geometries,  namely  in  petri 
dish  or  a  flow  cell,  in  full  3D  space  and  time. 

2.5.1  Boundary  conditions  for  biofilms  in  an  infinite  long  channel 

To  mimic  the  biofilm  development  in  a  culture  dish,  both  x  and  z  directions  are  assumed  periodic. 
At  boundaries  in  the  y  direction,  no-flux  boundary  conditions  are  imposed  for  the  biomass  and  functional 
components  and  no-slip  boundary  conditions  are  imposed  for  the  velocity  components,  that’s 


(ivs(f)s  -  Di4>sS7i)  ■  n\y=0Xy 

=  0, 

i  =  c,  d,  A , 

(2.26) 

^ 4*i  '  m\y=0,Ly 

=  0, 

i  =  61,  62,  63, 64,  p, 

(2.27) 

(' v<t>i  ~  A^V^-) '  n\y=o,Lv 

=  0, 

i  =  61,  62,  63,  64,  p. 

(2.28) 

We  also  impose  a  nutrient  feeding  condition  c\y=Ly  =  Co  in  place  of  the  zero-flux  condition  for  the  nutrient 
wherever  there  is  a  steady  supply  of  nutrients  through  the  boundary  at  y  =  Ly. 

2.5.2  Boundary  conditions  for  biofilms  in  a  finite  flow-cell 

For  the  case  of  biofilms  in  a  finitely  long  water  channel,  which  we  call  it  a  flow  cell,  periodic  boundary 
conditions  are  imposed  in  the  z  direction,  and  y  direction  is  bounded  by  solid  walls,  that’s  no  slip  boundary 
conditions  for  the  velocity  v  and  no  flux  boundary  conditions  for  the  biomass  and  functional  components 
are  imposed, 


(ivs(f)s  -  Di<t)s\/i )  •  n|j,=0,L„ 

=  0, 

i  =  c,  d,  A , 

(2.29) 

V 4>i  ■  Yl\y=0,Ly 

=  0, 

i  =  61,  62, 63,  64,  p, 

(2.30) 

1 

V 

< 

On  I  O 
^  |  *-+! 

3 

II 

0 

=  0, 

i  =  61,  62, 63,  64,  p. 

(2.31) 

The  channel  is  finite  in  the  x-direction.  So,  there  exist  the  inlet-outlet  boundary  conditions.  The  inlet  velocity 
at  x  =  0  is  given  by 

Vo  =  (poy(l  -  y),0,0) ,  (2.32) 

where  po  is  a  prescribed  pressure  gradient.  Suppose  that  the  solvent  has  already  reached  equilibrium  when 
flowing  out  of  the  cell  at  x  —  Lx,  we  impose  vx  =  0  at  the  out-let  boundary.  For  nutrient  c,  we  impose 
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the  feeding  boundary  is  at  x  =  0  with  a  prescribed  boundary  condition  c  =  co(y)  and  cx  =  0  is  assumed  at 
x  —  Lx.  For  the  biomass,  we  impose  no-flux  boundary  condition  in  the  x  direction,  which  only  applies  to 
the  situation  where  no  biomass  is  flown  in  and  out  of  the  flow-cell, 

-n\x=0}Lx  =0,  j  =  61, 62, 63,64, p.  (2.33) 

0(pi 


2.6  Nondimensionalization 


We  denote  to  as  the  reference  time  scale,  h  the  reference  length  scale,  and  co,  do  and  Aq  the  characteristic 
concentration  of  the  nutrient,  antimicrobial  agents  and  quorum  sensing  molecules,  respectively.  Then,  we 
nondimensionalize  the  variables  as  follows, 


t  „ 

t  —  — ,  x  = 
to 


tA  g-^o  ?  —  —  d  =  — 

Poh2  pohz  c0  do 


The  following  dimensionless  parameters  then  emerge, 


A  =  ^a,r 

to  ' 


7i  kTtQ  p 
poh 4  5  2 


!2kTtl 

poh 2 


.^.  =  “S.fep  =  £E.^  =  S£.P 


_  pph2 


_  poh 2 


Psto 


Ppto 


Pbto 


+ 


Pn 

PO 


Cbi  CbitOi  T bi  Cdz  Cai  ^z^O?  ^dz  ^dz^O?  ^  1,  2,  3. 


jZ  .  —  hbi  k  —  hdd  h  7  —  hd.  h  •  —  h  ■  —  ^EL  i  —  1  9  9 
I^bx  —  CQ  >  K'dz  —  d0  '  d  ~  d0'  ai  —  Ao'  P1  —  co  '  ~  ’  ’ 

fcp  =  rCpto,  kdp  =  ^,kcp  =  ^,Di  =  =  cO,  ^0,  -40. 


For  simplicity,  we  drop  the  symbol',  the  nondimentionalized  equations  in  the  hydrodynamic  model  for 
biofilms  are  summarized  as  follows: 


-5> 

+ 

c 

<1 

—  V  •  (( t>brb  +  rftp'hp  +  0s  Ts) 

-  Vp  -  TiV  •  (y4>nV(j)n 

V- V 

=  0, 

%i  +  V-(v06i) 

=  V  •  (A 05jV/ij)  +  j 

=  1,2, 3, 4, 

^  +  V-(v„^,) 

=  V  •  (A (j)pS7 Pp)  + 

^  +  V  •  (0svsc) 

=  V  •  (Dc</>sVc)  +  5c, 

%^  +  V-(v<M) 

=  V  •  (Da(/)S\7A)  +  ga, 

¥  +  V(^vsd) 

=  V  •  (Z)^0sVd)  +  gdi 

(2.34) 
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where  the  reactive  terms  are  given  by 
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X+l  ri4>bi  +  X)i=l  kj+d&bi  ~  J'4<^64> 

E?=1gfefc(i-si:)+W4 

-raJ4  +  (X^L,^  +  (l  -  -  *£+3^’ 

_  X+=l  k^r^^bi  -  X+=l  k^c^bii 

—  kd+d  (Si=l  rdi4)bi  T"  rdp(t>p )  —  rd^- 


rp(t)p  kfp+d^pi 


—  —i  a 


(2.35) 


3  Numerical  methods  and  GPU  implementation 


For  these  coupled  PDEs,  we  devise  a  numerical  method  to  solve  them.  In  each  time  step,  we  solve  the 
momentum  equation  first  using  the  extrapolated  data  for  the  biomass  and  functional  components.  Then  with 
the  updated  velocity,  we  solve  the  Cahn-Hilliard  equations  and  transport  equations  for  functional  compo¬ 
nents.  We  note  that,  in  the  following  notation,  any  variable  with  overline  72+1  represents  the  second-order 
extrapolation  from  the  n-lth  and  the  nth  to  the  n  +  1th  step. 

For  the  momentum  equation,  we  use  a  modified  Gauge-Uzawa  method  [?]  to  calculate  the  velocity. 
Recall  that  the  momentum  equation  is  given  by 

dv 

p(—  +  V  •  Vv)  =  V  •  {(f>bTb  +  (t>pTp  +  <f)sTs)  -  [Vp  +  TlV  •  (V0n  ®  V<pn)] 

By  adding  a  second  order  term  —  +-  V2v  on  both  sides,  where  Rea  is  a  numerically  chosen  Reynolds 
number,  we  rewrite  the  momentum  equation  into 

p{-^  +  v  •  Vv)  -  — v2v  =  -Vp  -  riV2<^nV(/>n  +  v  •  (< i>brb  +  (pPTP  +  (pSTS)  -  — v2v 


R 


ea 


Re 


The  modified  Gauge-Uzawa  method  is  given  by  three  steps. 
1.  Prediction: 


pn+l[MRR=^±Yl  *]  +  pn+l^n+1  .  V^n+l  +  d(V  •  (pn+1  vn+1))vn+1 , 

+  -+Vsn  +  Vpn  -  p^V2un+1  =  Rn+1  - 


u 


un+l 
n+1 1  _ 


n|j/=0,L„  =  0, 


=0  =  v0,  u: 


n+1 1 

x  \x 


=LX  =  0, 
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2.  Projection: 
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3.  Correction: 
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and  e  =  0.05  is  used  to  improve  the  stability  of  the  scheme.  The  numerically  chosen  Reynolds  number  R, 
is  computed  by 
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Where  Reb ,  RCp  and  i?6s  are  the  Reynolds  numbers  for  bacteria,  EPS  and  solvent,  respectively.  Here,  we 
note  that  5°  =  0  and  v1,  s1,  =  61, 62,  63,  64, p  and  c1,  d1,  A1  are  computed  by  a  first  order  scheme 

without  extrapolation. 

For  the  Cahn-Hilliard  equations,  we  use  a  semi-implicit  finite  difference  scheme.  Given  that 


4 
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we  can  easily  derive 
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Then  the  numerical  scheme  for  the  transport  equation  of  a  biomass  component  is  given  by 

''4  +,n+ 1  _l  An+ 1 
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(3.1) 


where  j  =  1,  2,  3, 4.  These  five  equations  are  solved  simultaneously  by  a  preconditioned  BI-CG-stab  rou¬ 
tine,  where  the  reactive  terms  are  extrapolated  to  the  n  +  1th  time-step.  The  reactive  terms  are  discretized 
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as  follows: 
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(3.2) 

At  each  time  step,  after  the  momentum  equation  and  the  transport  equations  for  biomass  components 
are  solved,  their  updated  data  would  be  used  to  solve  the  transport  equations  of  the  functional  molecules. 
These  equations  are  also  solved  by  a  series  of  second  order  semi-implicit  finite  difference  schemes,  which 
are  given  below 


30n+lcn+1— 4^nen+0n  — lcn— 1  +  ^+4  y^n+l^n+l)  =  y  .  (£>n+l$»+l y^+1)  +  gn+ 1} 
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+  vn+ 1  .  v(^4n++”+1)  =  V  •  (DJ+V^+1Van+1)  +  g2+1. 

where  the  discrete  schemes  for  reactive  terms  are  summarized  below 
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These  schemes  are  implemented  in  a  CPU-GPU  hybrid  environment  using  CUDA  to  achieve  paralleliza¬ 
tion  and  high  performance  computing.  The  implementation  has  been  well-tested  by  reducing  grid  sizes  in 
both  time  and  space,  and  nearly  second  order  convergence  rate  is  observed. 


4  Results  and  discussions 

The  hydrodynamic  model  of  biofilms  presented  in  section  2  is  fairly  general,  which  can  be  applied 
to  study  dynamic  problems  involving  biofilms  immersed  in  a  viscous  solvent.  In  this  section,  we  apply 
it  to  investigate  a  few  specific  issues  related  to  quorum  sensing  in  biofilm  formation,  which  include  the 
effect  of  the  threshold  for  quorum  sensing  induction  and  its  impact  to  biofilm  formation,  and  the  effect  of 
hydrodynamics  on  the  distribution  of  autoinducers  and  its  consequence  to  biofilm  development. 
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As  we  alluded  to  earlier,  the  model  contains  many  parameters,  some  of  which  need  to  be  calibrated  via 
experiments  while  others  are  either  obtained  from  the  literature  or  from  our  best  guesses  [?,  ?].  For  instance, 
the  wash-out-rates  in  the  model,  r\,  r<2,  r%,  r 4,  rp ,  ra  and  r d  are  normally  small;  so,  we  set  them  into  zero  in 
the  following  discussion.  For  the  monod  constants  (the  fc’s),  we  choose  a  common  value  3.5  x  10-4  as  used 
in  [?].  All  the  other  parameters  used  in  this  paper  are  summarized  in  table  ??,  unless  mentioned  otherwise. 


4.1  Reactive  kinetics  in  biofilms 


We  begin  with  the  study  on  reactive  kinetics  in  this  model  to  elucidate  bulk  reactive  kinetics  for  each 
component  on  two  important  mechanisms  during  biofilm  development.  By  neglecting  hydrodynamics,  dif¬ 
fusive  effects,  and  the  drug  effect  in  the  model,  as  well  as  assuming  nutrient  is  sufficiently  supplied,  we 
arrive  at  the  governing  equations  for  reactive  kinetics  in  biofilms, 
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(4.1) 


where  c—  1  is  used  in  the  following  discussion. 

It  is  observed  experimentally  that  biofilm  formation  can  be  affected  dramatically  by  the  threshold  for 
quorum  sensing  induction.  So,  in  the  first  case  study,  we  focus  on  investigating  the  solution  sensitivity  on 
the  threshold  parameter  r.  For  simplicity  and  without  loss  of  generality,  we  conduct  a  simulation  with  only 
quorum  sensing  bacteria  in  the  system,  where  we  simulate  biofilm  formation  with  a  few  selected  quorum 
sensing  threshold  values.  One  set  of  simulation  results  is  depicted  in  Figure  ??.  It  shows  that  the  lower  the 
threshold  value  is,  the  quicker  quorum  sensing  induction  can  take  place  so  that  quorum  sensing  up-regulated 
bacteria  can  begin  to  grow,  eventually  its  population  will  catch  up  or  even  overtake  that  of  the  quorum 
sensing  down-regulated  bacteria.  As  more  QS  up-regulated  bacteria  become  available,  the  amount  of  EPS 
production  increases,  leading  to  the  formation  of  robust  biofilms.  For  any  fixed  threshold  value,  a  steady 
state  exists  for  the  bacterial  populations  as  well  as  the  EPS  capacity  after  sometime  in  this  model.  When 
biofilm  development  reaches  a  steady  state,  the  ratio  of  the  concentration  of  the  up-regulated  bacteria  to  that 
of  the  down-regulated  bacteria  decreases  drastically  as  the  threshold  value  increases.  Besides,  the  transient 
exponential  growth  phase  for  bacteria  is  captured  and  the  transition  between  these  two  bacterial  phenotypes 
is  shown  at  the  moment  when  the  density  of  the  QS  molecules  reaches  the  threshold.  This  phenomenon  has 
been  observed  in  many  experiments  and  is  now  confirmed  by  this  simplified  bulk  kinetic  model.  We  remark 
that,  for  both  the  bacterial  growth  and  the  EPS  production,  we  implemented  logistic  growth  modes  in  the 
model  so  that  the  steady  state  is  limited  to  an  environmentally  sustainable  level. 

Secondly,  we  study  the  case  of  competition  between  QS  bacteria  and  non  QS  bacteria  in  a  biofilm  system 
with  respect  to  the  growth  rate  of  the  QS  down-regulated  bacteria  C&  1 .  There  are  many  parameters  in  this 
model,  our  extensive  parameter  studies  show  that  the  biofilm  dynamics  is  sensitive  to  the  growth  rate  C\}\  of 
the  QS  down-regulated  bacteria.  Suppose  that  we  have  10%  non-quorum  sensing  bacteria  initially  populated 
in  the  system.  We  simulate  biofilm  formation  with  respect  to  the  growth  rate  of  QS  down-regulated  bacteria. 
The  result  is  depicted  in  Figure  ??. 
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If  the  growth  rate  is  high,  the  total  bacterial  volume  fraction  is  not  affected  much  over  time  even  when  it 
reaches  a  steady  state,  and,  in  the  meantime,  non-quorum  sensing  bacteria  are  kept  in  a  low  volume  fraction. 
In  this  case,  the  biofilm  formation  with  an  exponential  growth  and  EPS  production  are  observed  in  transient. 
However,  when  the  growth  rate  is  low,  the  total  biomass  reaches  a  steady  state  at  a  slower  pace,  which 
is  primarily  controlled  by  non-quorum  sensing  bacteria  and  eventually  the  non  QS  regulated  bacteria  take 
control  over  the  whole  colony  should  the  initial  volume  fraction  of  the  non  QS  regulated  bacteria  be  high. 
The  result  shown  in  Figure  ??  is  with  respect  to  four  selected  values  of  growth  rate  C&i  in  6  days,  in  which 
the  volume  fraction  of  non-QS  bacteria  at  the  lowest  growth  rate  C&i  =  8  x  10-5  is  approximately  10  times 
higher  than  that  of  QS  bacteria  at  the  end  of  the  6th  day,  whereas  the  volume  fraction  of  non-QS  bacteria 
barely  changes  at  the  highest  growth  rate  of  C&i  =  8  x  10-4.  This  result  suggests  that  quorum  sensing 
in  biofilms  may  not  always  be  a  beneficial  mechanism  for  growth  when  there  is  a  large  amount  of  non-QS 
regulated  bacteria  present,  as  it  might  be  taken  advantage  of  by  other  phenotypes  or  even  enemies. 

A  complete  parameter  study  of  this  reactive  kinetics  model  perhaps  would  warrant  another  paper,  which 
will  not  be  pursued  further  in  this  paper.  We  next  turn  to  the  full  hydrodynamic  model  to  examine  two 
quorum  sensing  related  mechanisms  that  affect  the  heterogeneous  biofilm  formation. 

4.2  Biofilm  development  coordinated  by  quorum  sensing 

Biofilm  formation  is  a  process,  where  bacteria  aggregate  in  space  and  secrete  glue-like  exopopolysachar- 
rides  or  extra  cellular  polymeric  substances  (EPS)  to  form  a  bio-network  around  the  bacteria,  which  creates 
a  barrier  between  the  bacteria  and  the  surrounding  media.  It  has  been  observed  that  quorum  sensing  controls 
biofilm  formation  in  vibrio  cholerae  [?]  by  regulating  transcription  of  genes  involved  in  expopolysacharride 
production  (i.e.,  EPS  production).  Other  observations  on  quorum  sensing  regulated  EPS  production  can 
also  be  found  in  Pseudomonas  aeruginos  [?],  sinorhizobium  meliloti  [?],  as  well  as  in  some  other  bacterial 
systems  [?]. 

Since  there  has  not  been  a  consensus  on  what  is  the  role  of  quorum  sensing  during  biofilm  formation, 
we  will  focus  on  simulating  the  mechanism  of  quorum  sensing  regulated  EPS  production  using  the  model 
in  a  heterogeneous  biofilm  colony  and  studying  benefits  and  drawbacks  of  QS  regulation  during  biofilm 
formation  in  the  following. 

4.2.1  Biofilm  development  coordinated  by  quorum  sensing 

For  simplicity,  we  suppress  the  non-QS  regulated  bacteria  in  the  first  case  study,  i.e.,  we  consider  a 
model  biofilm  system  where  only  QS-regulated  bacteria  are  present  initially.  We  set  the  initial  value  of  3 
(the  volume  faction  of  non  QS  regulated  bacteria)  to  zero  and  assume  that  the  dead  bacteria  are  all  dissolved 
in  the  solvent  instantly.  This  way,  we  can  highlight  how  quorum  sensing  induction  takes  place  and  how 
the  biofilm  developmental  dynamics  can  be  coordinated  by  the  density  of  the  QS  molecules.  We  aim  our 
numerical  simulations  at  heterogeneous  biofilm  formation  in  an  aqueous  environment  where  laminar  flows 
can  be  present.  The  scenarios  include  a  culture  petri  dish  filled  with  water  or  a  long  water  channel  where 
nutrient  is  constantly  supplied  through  the  solution  surface. 

Figure  ??  depicts  a  result  from  the  simulation,  where  solutions  are  shown  at  t  —  0  and  t  =  300, 
respectively.  In  this  simulation,  some  QS  down-regulated  bacteria  ( QS~ )  are  assumed  to  be  attached  to 
the  substrate  surface  initially  and  the  bacterial  cell  reproduction  ensues  as  soon  as  the  simulation  begins. 
As  a  result,  the  QS  molecule  is  produced  and  its  concentration  keeps  growing.  Once  the  density  of  the 
quorum  sensing  molecules  exceeds  the  threshold  value  r,  QS~  start  to  bind  with  QS  molecules  and  undergo 
a  phenotypic  change  into  quorum  sensing  up-regulated  bacteria  QS+,  for  which  the  gene  expression  for 
reproduction  for  both  cells  and  EPS  is  switched  on.  The  biofilm  colony  starts  growing  exponentially  by 
reproducing  bacteria  as  well  as  secreting  EPS,  necessary  extracellular  RNA  and  enzymes.  With  more  QS 
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molecules  produced,  QS+  eventually  outgrow  QS~  and  populate  the  bulk  part  of  the  biofilm  colony.  In  the 
result,  we  observe  that  QS~  tend  to  be  seen  near  the  surface  of  the  biofilm  colony.  The  3D  result  for  volume 
fractions  of  total  bacteria,  QS  downregulated  and  QS  upregulated  bacteria,  EPS  and  the  total  biomass  are 
depicted  in  Figure  ??. 

The  model  predicts  that  QS -regulated  EPS  production  contributes  to  the  formation  of  heterogeneous 
biofilm  structures.  As  shown  in  Figure  ??(b),  the  distribution  of  live  bacteria  is  highly  heterogenous,  QS~ 
are  mainly  located  at  the  surface  of  the  biofilm,  where  they  can  access  nutrients  more  easily,  and  QS+  are 
located  more  in  the  bulk  of  the  biofilm  colony.  This  mushroom  shaped  biofilm  morphology  qualitatively 
agrees  with  the  published  result  on  the  morphological  pattern  in  Pseudomonas  aeruginosa  biofilms  [?]. 

In  order  to  show  details  of  this  process,  plots  of  three  slices  at  x  =  0.16,  y  =  0.04,  z  =  0.16  for  different 
biofilm  components  are  shown  in  Figure  ??  and  ??,  respectively.  In  Figure  ??(c),  where  the  average  velocity 
is  small,  the  distribution  of  the  QS  molecules  (AHFs)  is  dominated  mainly  by  spatial  diffusion.  Meanwhile, 
we  observe  that  QS  molecules  are  mainly  distributed  within  the  biofilm  colony  although  small  amount  can 
diffuse  into  the  solvent,  illustrated  by  a  three-slice  plot  in  Figure  ??(b).  The  biomass  flux  at  z  =  0.5  is 
plotted  in  Figure  ??(d).  It  shows  that  it  is  mainly  distributed  at  the  top  surface  between  the  biomass  the 
surrounding  solvent,  where  there  are  more  nutrients  available  at  the  surface  than  in  other  places  so  that  the 
biomass  tends  to  grow  the  fastest  in  these  places.  A  pressure  profile  is  also  plotted  in  a  slice,  which  shows 
that  a  high  pressure  is  observed  in  the  interior  while  a  low  one  is  present  near  the  surface.  This  numerical 
simulation  explains  the  growth  process  of  the  biofilm  due  to  QS  regulated  bacteria  and  QS-regulated  EPS 
production.  Of  course,  the  biomass  can’t  grow  without  a  sufficient  nutrient  supply.  But,  in  this  case  study, 
we  assume  its  role  is  secondary. 

4.2.2  Benefits  and  drawbacks  of  quorum  sensing  regulation 

Now  that  biofilm  formation  is  regulated  by  quorum  sensing,  a  straight  forward  question  is  if  this  quorum 
sensing  mechanism  is  always  beneficial?  Experimental  findings  seem  to  suggest  the  answer  is  ”it  depends.” 
As  the  production  of  QS  molecules  is  energetically  costly,  bacteria  probably  would  switch  on  the  quorum 
sensing  mechanism  to  start  producing  QS  molecules  whenever  they  are  in  the  ’’right  time”  and  ’’right  place.” 
In  this  subsection,  we  would  like  to  further  investigate  the  so  called  Tight  time’  and  Tight  place’  such  that 
the  QS  regulation  for  EPS  production  can  be  in  favor  of  the  biofilm  development. 

Through  numerical  simulations,  we  find  that  the  model  predicts  that  quorum  sensing  regulation  is  not 
beneficial  to  biofilm  formation  in  a  short  time  frame  when  competing  with  invading  strains  of  bacteria. 
In  other  words,  given  a  highly  reproductive  invasion  strain  (like  non-QS  regulated  bacteria),  QS-regulated 
bacteria  can  lose  the  competition  such  that  the  non-QS  regulated  bacteria  can  take  over  and  eventually 
populate  the  whole  biofilm  colony,  which  is  mentioned  in  the  case  study  of  bulk  reactive  kinetics. 

To  illustrate  the  point,  let’s  consider  a  biofilm  system  with  three  different  phenotypes  of  QS  regulated 
bacteria:  QS~ ,  QS+  and  non-QS".  In  this  system,  we  suppose  that  non-QS  regulated  bacteria  have  the 
same  reproductive  rate  as  that  of  QS~,  but  do  not  bind  QS  molecules  and  convert  themselves  into  QS+. 
We  conduct  a  numerical  simulation,  taking  into  account  all  these  three  phenotypes  to  study  the  detail  of 
how  the  interaction  between  these  bacteria  can  affect  biofilm  development  via  different  self-production 
and  EPS  production  rates.  In  this  simulation,  initially  some  QS~  and  non-QS  regulated  bacteria  with  the 
same  volume  fraction  are  assumed  to  be  attached  to  the  substrate  surface.  We  assume  C&i  =  C&3  and  the 
growth  rate  for  QS  up-regulated  bacteria  C&2  is  smaller.  A  3D  view  of  the  biofilm  growth  process  predicted 
by  our  3D  simulations  is  displayed  in  Figure  ??.  As  shown  in  Figure  ??(b),  non-QS  regulated  bacteria 
outgrow  the  QS  regulated  bacteria  and  take  over  the  whole  biofilm  colony  eventually.  This  is  because  the 
QS  regulated  bacteria  undergo  mutual  conversion  which  slows  down  the  reproduction  of  QS+  bacteria  as 
well  as  QS~ .  The  total  volume  fraction  for  each  type  of  bacteria  is  plotted  in  ??(c),  where  non-QS  bacteria 
grow  exponentially,  while  QS  regulated  bacteria  grows  only  slightly.  Since  non-QS  regulated  bacteria  have 
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a  higher  effective  reproduction  rate  than  those  of  QS  regulated  ones  combined  effectively,  non-QS  regulated 
bacteria  eventually  grow  faster  to  gain  more  access  to  nutrient.  Hence,  the  density  of  non-QS  regulated 
bacteria  can  eventually  take  over  that  of  QS  regulated  bacteria  in  a  majority  part  of  the  biofilm  colony  due 
exclusively  to  the  quorum  sensing  effect. 

In  addition,  our  simulations  also  show  that  quorum  sensing  regulation  is  not  beneficial  to  biofilm  devel¬ 
opment  in  a  short  time,  especially,  in  terms  of  gaining  biomass.  To  expose  the  feature,  we  conduct  another 
numerical  simulation,  where  the  quorum  sensing  mechanism  is  switched  off  by  setting  7*21  =  r\2  —  0  while 
the  other  model  parameters  are  retained  as  in  the  previous  numerical  experiment.  It  is  shown  in  Figure  ?? 
that  the  volume  fraction  of  the  total  live  bacteria  without  quorum-sensing  regulation  is  higher  than  that  with 
quorum  sensing  regulation  switched  on.  However,  the  total  biomass  in  the  system  without  QS  regulation, 
including  EPS  and  bacteria,  is  much  less  than  that  with  QS  regulation.  This  indicates  that  QS  regulation 
is  not  a  way  to  reproduce  bacteria,  but  rather  a  way  to  increase  the  amount  of  EPS  so  as  to  expand  biofilm 
colony,  especially,  in  a  short  time  frame.  The  result  agrees  qualitatively  with  the  experimental  observation 
reported  in  [?]. 

In  the  long  term  however,  our  numerical  studies  show  that  QS  regulation  is  of  benefit  to  stabilizing 
biofilm  colony  and  settling  on  the  host,  as  well  as  protecting  bacteria  from  environmental  stress.  Given  the 
protective  nature  of  EPS,  the  biofilm  colony  can  become  more  stable  and  resistant  to  external  stress  while 
more  EPS  is  produced.  EPS  can  protect  bacteria  from  being  attacked  by  antimicrobial  agents  or  immune 
systems  of  the  host.  To  illustrate  the  beneficial  effect  of  quorum  sensing,  we  conduct  a  numerical  experiment 
with  the  drug/disinfectant  effect  included.  The  result  is  shown  in  Figure  ??.  In  this  simulation  (see  Figure 
??(f)),  non-QS  regulated  bacteria  are  killed  dramatically  while  QS  regulated  bacteria  are  killed  in  a  much 
slower  pace.  Our  model  and  simulation  confirm  that  EPS  acts  as  a  barrier  to  protect  bacteria  from  being 
attacked.  In  Figure  ??(a)  and  Figure  ??(b),  we  notice  that  bacteria  staying  within  the  environment  filled 
with  EPS  lives  longer  under  antimicrobial  drug  treatment.  In  Figure  ??(e),  we  see  that  drug  penetrates  less 
in  the  place  where  there  exists  more  EPS  than  in  the  place  with  less  EPS.  This  numerical  result  agrees  qual¬ 
itatively  with  the  claim  made  in  [?],  that  biofilms  can  benefit  from  quorum  sensing  induced  EPS  production 
if  bacteria  cells  have  the  objective  of  acquiring  a  thick,  protective  layer  of  EPS.  This  may  also  explain  why 
antimicrobial  drugs  are  more  effective  for  young  biofilms;  when  biofilms  have  grown  for  weeks,  the  effect 
of  antimicrobial  drugs  drops  significantly.  When  the  biofilm  is  young,  there  exists  less  EPS  surrounding  it; 
when  the  biofilm  is  ’’old”,  however,  the  QS  induced  EPS  has  already  accumulated  around  bacterial  cells  so 
that  they  can  prevent  drugs  from  penetrating  deeper  into  the  biofilm.  This  agrees  qualitatively  well  with  the 
observation  reported  in  [?]  that  biofilms  at  an  earlier  age  is  more  susceptible  to  antimicrobial  treatment  than 
those  in  an  older  age. 

We  have  used  two  numerical  experiments  to  elucidate  the  role  of  quorum  sensing  in  the  development  of 
biofilms.  Benefits  or  detriments  of  quorum  sensing  for  biofilms  depend  on  many  factors,  for  instance,  the 
age  of  biofilms,  if  antimicrobial  treatment  is  involved,  growing  the  biofilm  or  eradicating  the  biofilm,  etc. 
Quorum  sensing  regulation  tends  to  slow  down  the  production  of  bacteria  while  in  the  meantime  promote  the 
production  of  EPS,  which  is  beneficial  to  the  development  of  a  robust  biofilm  in  a  long  run  but  detrimental 
to  the  biofilm  formation  in  a  short  time.  Nevertheless,  quorum  sensing  is  an  important  mechanism  for 
biofilm  development.  An  ultimate  determination  of  its  role  can  only  be  rendered  in  connection  with  the 
specific  situation  where  biofilm’s  issues  are  concerned.  Next,  we  turn  to  discussing  the  interaction  between 
hydrodynamics  and  quorum  sensing  during  the  development  of  biofilms  in  an  aqueous  environment. 
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4.3  Interaction  between  biofilm  hydrodynamics  and  quorum  sensing 

4.3.1  Quorum  sensing  in  an  aqueous  environment 

Many  biofilms  grow  in  an  aqueous  environment  ranging  from  laminar  flows,  such  as  river  bank  or  pipes, 
to  quiescent  aqueous  environment.  It  is  observed  that  a  flow  environment  can  impact  on  the  QS  induction 
during  biofilm  formation  [?,  ?].  In  return,  the  hydrodynamically  altered  QS  induction  can  affect  biofilm 
structures  and  functions.  However,  there  is  little  known  on  the  detail  of  how  hydrodynamic  factors  actually 
affect  QS  induction.  Here,  we  use  the  hydrodynamic  model  along  with  the  developed  numerical  simulating 
tool  to  investigate  some  of  the  effects  in  details.  Again  for  simplicity,  we  consider  the  case  there  are  only 
QS  regulated  bacteria  in  the  biofilm  system,  i.e.  QS~  bacteria  and  QS+  bacteria,  as  we  are  interested  in 
how  quorum  sensing  induction  is  affected  by  hydrodynamic  factors.  The  boundary  conditions  used  in  the 
simulation  are  set  up  to  mimic  biofilm  development  in  a  flow  cell,  where  the  fluid  (solvent)  flows  through 
continuously. 

Firstly,  we  conduct  a  numerical  simulation  of  the  biofilm  system  subject  to  a  varying  inlet  fluid  velocity 
and  focus  on  biofilm  formation  under  hydrodynamic  shear.  The  result  is  shown  in  Figure  ??.  Initially,  some 
QS~  bacteria  are  attached  at  the  center  of  the  flow  cell.  Solvent  containing  nutrient  keeps  flowing  through 
the  flow-cell  while  the  bacteria  grow  into  a  biofilm  colony.  As  shown  in  Figure  ??(c),  a  higher  hydrodynamic 
shear  can  actually  dilute  the  concentration  of  QS  molecules,  which  leads  to  a  delay  in  QS  induction.  With  the 
diluted  QS  molecular  concentration,  the  biofilm  in  the  flow  cell  with  a  higher  inlet  velocity  has  a  less  volume 
fraction  of  QS+  bacteria.  Note  however  that  a  biofilm  with  a  higher  inlet  velocity  results  in  a  higher  volume 
fraction  of  the  total  bacteria.  This  is  because  the  nutrient  supply  at  a  higher  inlet  velocity  can  bring  in  more 
nutrient  and  thus  facilitates  the  growth  of  all  bacteria.  A  robust  biofilm  colony  needs  not  only  the  growth  of 
bacteria,  but  also  the  build-up  of  EPS.  Since  the  EPS  growth  by  QS+  is  more  than  that  by  QS~ ,  a  reduced 
QS+  production  leads  to  the  reduced  EPS  production  in  the  biofilm.  Even  though  the  total  bacteria  is  not 
reduced,  the  EPS  production  is  reduced  by  the  higher  in-let  velocity.  2D  slices  (at  z  =  0.5)  of  QS  molecules 
and  QS+  bacteria  at  time  t  —  25  are  shown  in  Figure  ??.  This  numerical  experiment  demonstrates  that  the 
maximum  concentration  of  Qs  molecules  and  the  volume  fraction  of  QS+  are  higher  in  the  flow  cell  when 
subject  to  a  slower  inlet  speed  and  the  scenario  is  reversed  when  the  inlet  speed  is  higher. 

This  study  suggests  that  hydrodynamics  can  have  a  strong  influence  on  quorum  sensing.  A  higher  inlet 
flow  dilutes  the  concentration  of  the  QS  signaling  molecules,  which  then  requires  a  higher  cell  density  and 
delayed  time  for  quorum  sensing  induction  to  take  place.  This  means  that  QS  may  not  be  fully  operational 
at  high  flow  rates  and  thus  indicates  that  the  role  of  QS  in  biofilm  formation  can  be  compromised  in  a  flow 
environment. 

We  note  that  how  quickly  QS  molecules  are  produced,  how  quickly  the  QS  molecules  diffuse  away 
from  the  biofilms  and  the  critical  concentration  of  QS  molecules  are  three  main  factors  affecting  biofilm 
formation.  When  taking  into  account  of  the  effect  of  nutrient,  especially,  when  the  growth  due  to  nutrient 
dominates,  the  diffusion  rate  and  consumption  rate  of  nutrient  can  be  essential  for  the  biofilm  growth  as 
well.  We  hope  the  theoretical  result  can  be  quantitatively  confirmed  experimentally  in  the  near  future. 

4.3.2  Effect  of  nutrient  to  the  interaction  between  hydrodynamics  and  quorum  sensing 

Compared  with  a  quiescent  aqueous  environment,  the  big  difference  of  biofilms  developed  in  a  flow  cell 
is  their  structural  heterogeneity  in  the  flow  field,  due  to  the  hydrodynamic  stress  and  inhomogeneous  nutrient 
distribution.  Mass  transfer  in  the  fluid  mixture  can  be  affected  by  hydrodynamics  of  the  bulk  fluid  and  the 
morphology  of  the  biofilm  colony.  Under  a  flow  environment,  one  natural  question  to  ask  is  which  part  of 
biofilms,  especially  bacteria  living  upstream  or  downstream,  is  benefited  from  quorum  sensing  regulation 
more.  For  this,  we  conduct  a  numerical  simulation  with  a  distribution  of  biofilm  colonies  in  the  flow  channel. 
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Our  numerical  simulations  suggest  that  biofilms  upstream  can  induce  downstream  biofilms  to  quorum- 
sense  under  the  flow  if  they  are  close  enough  spatially,  which  is  in  agreement  with  the  result  in  [?].  In 
Figure  ??,  we  demonstrate  that  quorum  sensing  takes  place  at  downstream  firstly,  that’s  quorum  sensing 
induction  takes  place  at  the  last  four  colonies  instead  of  the  first  one  located  upstream.  It  also  shows  that  the 
downstream  region  tends  to  have  a  higher  concentration  of  QS  molecules  than  at  the  location  upstream  in 
Figure  ??(g),  which  we  believe  is  an  amplification  of  the  convection  and  diffusion  effect  of  Qs  molecules  . 

However,  we  note  that  this  may  not  be  always  true,  especially  when  nutrient  supply  is  strong.  A  compar¬ 
ison  with  different  nutrient  supply  rates  is  shown  in  Figure  ??.  We  observe  that,  at  an  early  stage,  bacteria 
downstream  can  be  benefitted  since  QS  molecules  are  more  diluted  upstream  (see  Figure  ??(c)).  QS  in¬ 
duction  is  more  effective  downstream  since  there  exist  more  QS+  bacteria  there  (see  Figure  ??(b)).  But 
when  time  passes  by  and  the  nutrient  supply  is  kept  at  a  sufficiently  high  level,  it  seems  that  quorum  sensing 
induction  benefits  from  a  higher  flow  velocity  since  it  sustains  a  sufficient  nutrient  supply  that  facilitates 
the  growth  of  all  bacteria  as  well  the  production  of  QS  molecules.  As  a  result,  the  overall  production  of 
QS  molecules  upstream  overcomes  the  dilution  of  QS  molecules  such  that  both  the  concentration  of  QS 
molecules  and  the  volume  fraction  of  bacteria  are  higher  upstream.  In  Figure  ??(e),  we  can  see  there  are 
more  QS+  bacteria  upstream  than  downstream  in  this  case.  Figure  ??(h)  depicts  the  EPS  production  in  the 
flow  cell,  which  can  also  be  tied  to  the  excessive  accessibility  to  nutrient. 

4.3.3  Biofilm  morphology  under  a  strong  shear  flow 

In  previous  discussions,  we  focus  on  interactions  between  hydrodynamics  and  quorum  sensing  induction 
under  the  assumption  that  the  biofilm  grows  under  weak  flow  fields.  When  the  flow  is  strong,  however, 
biofilms  can  undergo  significant  morphological  changes  and  most  of  the  biofilms  can  perhaps  be  flushed 
out  of  the  flow  cell.  To  better  illustrate  the  hydrodynamic  impact  on  biofilm  structures  and  QS  induction  in 
strong  flows,  we  conduct  a  numerical  simulation  where  a  grown  biofilm  colony  is  subject  to  a  string  shear 
flow  in  a  flow  cell.  Here,  we  put  two  biofilm  colonies  in  a  4  x  1  x  1  domain,  representing  a  flow  cell  with  inlet 
and  outlet  boundaries  in  the  x  direction.  The  initial  biofilm  profile  is  shown  in  Figure  ??(a).  From  ??(b), 
we  notice  that  the  biofilm  undergoes  a  significant  morphological  change  under  the  shear.  The  colonies  are 
tilted  under  the  shear.  The  one  near  the  inlet  boundary  is  about  to  detach  from  the  bottom  in  Figure  ??(d). 
In  addition,  a  fast  flow  can  dilute  QS  molecules  as  they  are  readily  convected  out  of  the  flow  cell  shown  in 
Figure  ??(c)  and  (e). 

The  advantage  of  our  full  3D  hydrodynamic  model,  compared  with  those  in  lower  dimension,  as  well 
as  reaction-diffusion  models  is  its  capability  to  resolve  the  flow  detail  and  hydrodynamically  as  well  as 
Theologically  significant  quantities  such  as  the  normal  stress  differences  and  shear  stress  components.  Here, 
these  stress  components  related  to  the  simulation  above  are  depicted  in  Figure  ??.  In  Figure  ??,  we  only 
show  the  sign  of  </>nr,  the  effective  stress  for  the  biomass  component,  where  we  use  red  to  represent  negative 
and  blue  positive  values  of  the  stress.  It  is  seen  through  Figure  ??(e)  that  the  first  normal  stress  difference 
is  positive  at  the  top  of  the  biofilm  region  while  negative  around  the  neck  facing  the  inlet  boundary,  which 
means  that  the  normal  stress  rxx  is  larger  than  ryy  at  the  top  of  the  biofilm  mechanically  providing  a  drag 
force  along  the  x  direction.  This  is  confirmed  by  Figure  ??(c),  where  the  top  of  the  biofilm  is  stretched 
towards  the  x  direction.  In  addition,  through  Figure  ??(D),  the  second  normal  stress  difference  is  positive  at 
the  neck  facing  the  inlet  boundary,  where  ryy  is  larger  than  rzz  resulting  in  a  stretching  of  the  biofilm  colony 
toward  the  y  direction.  This  is  supported  by  the  fact  that  the  neck  becomes  thinner  due  to  this  stretch,  shown 
in  Figure  ??(c).  In  addition,  the  other  three  components  of  shear  stress  tensor  are  also  depicted  in  Figure 
??(f-g).  On  a  y  plane,  i.e.  the  plane  y  =  a  (0  <  a  <  1),  we  see  the  shear  stress  rxy  is  positive,  which  forces 
the  biofilm  towards  the  x  direction  (shown  in  the  comparison  of  Figure  ??(b)  and  (c).)  The  effect  of  rxz  and 
ryz  can  also  be  observed  on  a  z-plane,  i.e.  z  =  a  (0  <  a  <  1),  shown  in  Figure  ??(g),  where  rxz  provides 
a  shear  stress  towards  the  x  direction  on  the  top  of  the  biofilm  mushroom  shaped  colony  and  ryz  provides  a 
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shear  stress  towards  the  y  direction  on  the  top  of  the  mushroom. 


5  Conclusion 

In  this  paper,  a  hydrodynamic  phase  field  model  for  biofilms  is  developed  to  study  biofilm  development 
coordinated  by  quorum  sensing,  extending  our  previous  biofilm  models  for  binary  mixture  [?].  By  identify¬ 
ing  appropriate  bacterial  types  based  on  their  function  in  quorum  sensing  in  the  biofilm  and  proposing  the 
corresponding  reactive  kinetics,  we  are  able  to  describe  several  bulk  phenomena  related  to  quorum  sensing 
induction  using  the  new  model.  We  then  develop  a  numerical  scheme  to  solve  the  governing  system  of 
equations  in  the  model  and  implement  it  in  3D  space  and  time.  Then  we  are  able  to  simulate  and  predict 
heterogeneous  spacial-temporal  structures  of  a  biofilm  in  terms  of  its  volume  fractions  for  each  bacterial 
component,  EPS  and  solvent.  The  simulation  tool  is  put  in  use  to  study  a  set  of  selected  mechanisms  in 
quorum  sensing  regulated  biofilm  development. 

Our  numerical  simulations  indicate  that  QS-regulated  EPS  production  acts  as  a  means  of  competition 
for  both  space  and  nutrient  since  higher  EPS  production  rate  can  help  to  expand  biofilm  colonies  spatially 
without  reaching  a  higher  local  density.  In  addition,  numerical  simulations  confirm  the  hypothesis  that  QS 
regulated  EPS  production  acts  as  a  means  for  promoting  biofilm  development  as  well  as  protecting  biofilms 
from  the  environment  stress,  such  as  antimicrobial  treatment  and  hydrodynamical  stress  etc.  However, 
quorum  sensing  regulation  is  not  universally  beneficial  to  biofilm  development  in  all  situations  evidenced 
by  what  we  have  shown  that  an  invasion  strain  of  bacteria  with  a  high  production  rate  can  out-grow  QS 
regulated  bacteria,  as  QS  induction  can  retard  the  proliferation  of  bacteria,  thereby  the  development  of 
biofilms,  in  a  short  time  window. 

By  simulating  biofilms  in  a  flow  cell,  hydrodynamic  effects  on  quorum  sensing  are  investigated  in  de¬ 
tails.  The  model  can  capture  the  competition  between  the  dynamic  diluting  effect  and  the  nutrient  supply 
for  QS  induction  well.  Our  simulations  predict  that  downstream  bacteria  can  be  benefited  from  higher  in 
flow  speed  to  initiate  QS  induction  earlier.  However,  the  model  also  shows  that  this  only  takes  place  in  a 
short  time  window,  where  the  nutrient  supply  is  weak.  When  the  nutrient  supply  is  strong,  it  is  a  competi¬ 
tion  between  the  hydrodynamic  dilution  of  QS  molecules  and  the  nutrient  supply  for  the  production  of  QS 
molecules  that  ultimately  determines  which  part  of  biofilms  in  a  flow  cell  would  benefit  more. 

This  3D  hydrodynamic  model  provides  a  powerful  predicative  tool  to  allow  us  to  investigate  quorum 
sensing  regulated  biofilm  development  and  some  factors  that  correlate  with  quorum  sensing  and  their  inter¬ 
action  with  hydrodynamics.  It  is  a  robust,  hands-on  tool  for  analyzing  the  mechanism  of  quorum  sensing  in 
biofilm  dynamics.  In  addition  to  the  phenomena  that  we  numerically  confirmed,  we  hope  results  predicted 
here  can  be  verified  experimentally  in  the  near  future.  In  our  subsequent  research,  other  quorum  sensing 
regulated  factors  such  as  virulence,  resistance  to  antimicrobial  agents  will  be  investigated. 
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Table  1:  Parameters  and  their  Values.  Parameters  with  references  for  their  order-of-magnitude 
are  marked,  otherwise  they  are  fitted  by  our  previous  work  and  experience. 


Symbol 

Description 

value 

Unit 

Source 

h 

Characteristic  length  scale 

1  x  10-3 

m 

[?] 

to 

Characteristic  time  scale 

1000 

s 

[?] 

Lx  5  Ly  ,  Lz 

size  of  computational  domain 

1  -  2  x  lO"3 

m 

[?] 

co,  do,  ao 

Characteristic  Nutrient  concentration 

8.24  x  10“3 

kg  m~3 

[?] 

T 

Absolute  Temperature 

303 

K 

[?] 

ks 

Boltzmann  constant 

1.38  x  10“23 

m2kgs~2  K~ 

1  [?] 

71 

Distortional  energy  coefficient 

8  x  106 

m~ 1 

[?] 

72 

Mixing  free  energy  coefficient 

3  x  1017 

m~3 

[?] 

X 

Flory  Huggins  parameter 

0.55 

[?] 

A 

Mobility  parameter 

1  x  10“9 

kg~1m3s 

N 

Generalized  polymerization 

1  x  103 

[?] 

Pn 

Network  density 

1  x  103 

kg  m~6 

[?] 

Ps 

Solvent  density 

1  x  103 

kg  m~3 

[?] 

rip,  Pb 

Dynamic  viscosity  of  biomass 

10 

kg  m~1s~ 1 

Ps 

Dynamic  viscosity  of  solvent 

1.002  x  10“3 

kg  m~1s~ 1 

[?] 

Dco 

Nutrient  diffusion  coefficient 

2.3  x  10“a 

m~2s~i 

[?] 

DaO 

AHL  diffusion  coefficient 

2.3  x  lO"10 

m~2s~ 1 

[?] 

Ddo 

Antimicrobial  diffusion  coefficient 

2.3  x  10“9 

m~2s~1 

[?] 

chl 

bacteria  growth  rate  for  QS~ 

4.0  x  10-4 

s~L 

[?] 

Cb2 

bacteria  growth  rate  for  QS+ 

1.0  x  10“4 

s_1 

[?] 

cb  3 

bacteria  growth  rate  for  non-QS 

4.0  x  10“4 

S_1 

[?] 

Cdl,  Cd2,  Cd3 

death  rate  due  to  drug 

4.0  x  lO”2 

s_1 

[?] 

r  12,  r*2i 

transfer  rate  between  foi  and  4>b2 

4  x  10"3 

s_1 

Cpl  ,  Cp3 

Growth  rate  of  <frp  due  to  (f)bi  and  4>b 3 

1  x  10“Y 

CP2 

Growth  rate  of  cf)p  due  to  cf)b2 

4  x  10“4 

5-1 

r’al,  ra2.  Cal,  Ca2 

growth  rate  of  AHL 

4  x  10“5 

s_i 

r di,  rdp,  rbi,  rpi,  i  =  1,2,3 

drug  and  nutrient  consumption  rate 

4  x  10-* 

s_i 

Dpr 

Hinson  model  constant 

0.02 

[?] 

fib, max 

carrying  capacity  for  bacteria 

0.15 

(ftp,  max 

carrying  capacity  for  EPS 

0.2 

T 

threshold  of  quorum  sensing 

0.5 

P 

bacteria  growth  rate  ratio 

1.0 
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tt 
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Time  (days) 

(a)  <j>b\ 


(b)  d>b2 


Figure  2:  Biofilm  development  with  only  QS  bacteria  in  the  biofilm  system.  This  simula¬ 
tion  demonstrates  the  effect  of  QS  threshold  value  r  for  quorum  sensing  induction  on  biofilm 
development.  The  QS  threshold  values  used  are  r  =  0.3, 0.5, 0.7, 0.9,  respectively,  and 
ri2  —  r2i  —  10~4,/?  =  0.  The  volume  fraction  of  QS  down-regulated  bacteria,  QS  up- 
regulated  bacteria,  EPS  and  total  bacteria,  respectively,  and  AHL  concentration  are  shown  as 
functions  of  time  in  the  subfigures,  (a)  The  volume  fraction  of  QS  down-regulated  bacteria, 
(b)  The  volume  fraction  of  QS  up-regulated  bacteria,  (c)  The  volume  fraction  of  EPS.  (d)  The 
volume  fraction  of  the  total  bacteria,  (e)  The  concentration  of  autoinducers/QS  molecules: 
AHL.  The  simulation  shows  that  the  smaller  the  threshold  value  r  is,  the  sooner  QS  induction 
takes  place.  The  sooner  QS  induction  takes  place,  the  faster  QS  up-regulated  bacteria  can  grow, 
leading  to  steady  states  in  EPS  production  as  well  as  the  total  biomass  sooner.  The  production 
of  QS  molecules  decreases  as  the  threshold  value  increases. 
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Figure  3:  Biofilm  development  with  both  QS  and  non-QS  bacteria  in  the  biofilm  system.  The 
numerical  simulation  begins  with  90%  QS  bacteria  and  10%  non-QS  bacteria  among  all  bac¬ 
teria  whose  volume  fraction  is  2%.  Four  values  of  the  growth  rate  for  down-regulated  QS 
bacteria  C&i  are  tested:  C& i  =  8  x  10-5,  2  x  10-4,  4  x  10-4  and  8  x  10-4  with  (3  =  0. 
The  corresponding  volume  fraction  for  each  component  is  shown  in  the  figure,  (a)  The  volume 
fraction  of  QS-down  regulated  bacteria,  where  a  rapid  growth  occurs  initially,  it  is  followed 
by  a  mild  decay  in  a  plateau,  then  a  rapid  decay,  and  finally  the  volume  fraction  settles  down 
to  a  near  steady  state;  (b)  the  volume  fraction  of  QS-up  regulated  bacteria,  which  increases 
monotonically;  (c)  the  volume  fraction  of  non-QS  bacteria,  which  increases  rapidly  to  a  steady 
state  value;  (d)  the  volume  fraction  of  the  total  biomass  (sum  of  the  volume  fractions  of  bac¬ 
teria  and  EPS);  (e)  the  concentration  of  QS  molecules,  AHL;  (f)the  volume  fraction  of  EPS. 
The  simulation  shows  that  the  higher  C&i  is,  the  faster  QS  bacteria  grow,  which  suppresses  the 
growth  of  non-QS  bacteria  and  enables  the  total  biomass  to  reach  a  steady  state  sooner.  At  a 
small  value  of  C&i,  the  growth  of  QS  bacteria  is  slow  which  allows  the  reproduction  of  non-QS 
regulated  bacteria  to  outgrow  QS  regulated  bacteria  eventually.  In  this  case,  it  also  delays  the 
biomass  to  reach  a  steady  state  so  that  the  entire  process  of  forming  a  robust  biofilm  is  slowed 
down. 
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Figure  4:  Biofilm  formation  with  EPS  production  coordinated  by  quorum  sensing  in  the  biofilm 
system  where  non  QS  bacteria  are  neglected.  This  simulation  shows  several  clusters  of  QS~ 
in  a  developing  heterogenous  biofilm  colony,  where  EPS  production  during  biofilm  formation 
is  coordinated  by  quorum  sensing.  The  initial  profile  is  shown  in  (a).  The  characteristic  time 
scale  used  is  to  =  103  seconds  and  the  simulation  is  conducted  from  t  =  0  to  t  =  300,  which 
represents  roughly  3.4  days  in  real  time.  (b).  A  time  t  =  300,  the  bacterial  distribution  in  terms 
of  volume  fractions  is  shown,  where  the  red  represents  QS~  and  the  blue  represents  QS+\  they 
are  mixed.  In  the  result,  more  QS+  are  observed  in  the  bulk  of  the  biofilm  while  more  QS~ 
show  up  near  the  biofilm  surface.  In  (c)-(f),  three  2D  slices  of  QS~ ,  QS+,  EPS  and  the  total 
biomass  at  t  =  300  are  depicted,  respectively,  from  which  more  interior  details  of  the  biofilm 
colony  are  on  display. 
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(a)  Nutrient  at  t=300 


(b)  AHL  at  t=300 


(c)  v  at  t=300 


(d)  0nvn  at  t=300  (e)  Pressure  at  t  —  300 

Figure  5:  Velocity,  pressure  and  functional  molecules  from  the  simulation  shown  in  Figure  ??. 
(a)-(b)  2D  slices  for  the  distribution  of  the  nutrient  and  QS  molecules,  AHL;  (c)-(e)  2D  slice 
at  z  =  0.5  for  the  average  velocity  v  and  the  biomass  flux  as  well  as  the  hydrostatic  pressure 
at  t  =  300,  respectively.  Low  pressures  occur  in  the  interfacial  region  where  the  biofilm  tends 
to  grow,  (d)  depicts  explicitly  the  biomass  flux  given  by  This  shows  that  QS  molecules, 

AHL,  are  mainly  concentrated  in  the  biofilm  colony,  but,  small  amount  is  diffused  into  the 
surrounding  environment. 
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(a)  Initial  Condition 


(b)  Bacteria  at  t=252  (c)  Volume  fractions  for  each  bacterial 

type 


(d)  063  at  t=252 


(e)  06i  and  062  at  t=252 


(f)  063  at  t=252 


(g)  06i  and  062  at  t=252 


(h)  Pressure  at  t=252 


(i)  Biomass  flux  0nvn  att=252 


Figure  6:  Biofilm  development  with  both  QS  bacteria  and  Non-QS  bacteria  in  the  biofilm  sys¬ 
tem.  Randomly  distributed  biofilm  colonies  of  multiple  bacterial  types  are  given  initially.  Here 
the  characteristic  time  is  set  at  to  =  103  seconds,  (a)  Initial  bacterial  volume  fraction  profiles; 
(b)  bacterial  volume  fractions  at  time  t  =  252,  where  purple  represents  non-QS  bacteria  and 
red  represents  QS-bacteria;  (c)  the  volume  fraction  for  each  bacterial  phenotype  in  the  entire 
biofilm  colony,  where  non  QS  regulated  bacteria  outgrow  the  QS  regulated  bacteria;  (d)the 
volume  fraction  of  non  QS  bacteria;  (e)  the  volume  fraction  of  QS  bacteria  at  t  =  252;  (f)-(i) 
2D  slices  at  z  =  0.5  for  non-QS,  QS  bacteria,  hydrostatic  pressure,  and  biomass  flux  </>nvn 
at  t  =  252,  respectively.  The  non  QS  regulated  bacteria  outgrow  the  QS  regulated  bacteria  in 
most  part  of  the  biofilm  except  in  regions  near  the  bottom  of  the  domain.  The  pressure  is  again 
small  near  the  interface  and  large  in  the  interior  of  the  biofilm. 
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Figure  7:  Comparison  between  bacterial  volume  fractions  during  biofilm  development  for 
biofilms  with  and  without  quorum  sensing.  Here  we  set  r\2  =  7*21  =  0  for  simplicity,  (a). 
Volume  fractions  of  total  biomass,  (b).  Volume  fractions  of  QS  down-regulated  bacteria,  (c). 
Volume  fractions  of  QS  up-regulated  bacteria,  (d).  Volume  fractions  of  total  live  bacteria.  The 
initial  conditions  are  the  same  as  the  one  used  in  Figure  ??.  With  quorum  sensing,  the  total 
biomass  is  higher  than  without  quorum  sensing,  but  the  total  live  bacteria  is  not  necessary  be 
higher. 
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(c)  Dead  Bacteria  at  t=94 


(d)  Three  Slices  of  4>bA  at  t=94 


x  10 


(e)  Three  Slices  of  Drug  at  t=94 


Figure  8:  Biofilm  development  subject  to  both  quorum  sensing  and  antibiotic  (drug)  treatment. 
The  initial  biofilm  profile  is  a  grown  biofilm  at  t  —  252  from  Figure  ??.  Antibiotic  agents  are 
dosed  from  the  top  boundary.  Here  the  characteristic  time  is  to  =  10  seconds,  (a)-(c)  shows 
the  volume  fraction  distribution  of  the  live  bacteria,  EPS,  dead  bacteria  at  t  —  94,  respectively; 
(d)-(e)  shows  three  slices  of  dead  bacteria  and  antibiotics  (drug)  at  t  =  94,  respectively;  (f) 
shows  total  volume  fractions  of  QS  bacteria  and  non-QS  bacteria  from  t  —  0  to  t  =  600, 
respectively.  After  a  long  period  of  antibiotic  treatment,  the  bacterial  population  is  reduced  to 
a  negligible  level  in  this  simulation. 
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(c)  Total  volumes  of  each  components 

Figure  9:  Biofilm  development  under  different  inlet  flow  velocities.  This  shows  a  comparative 
study  on  how  some  QS~  bacteria  develop  into  a  biofilm  colony  under  different  inlet  velocities 
or  equivalently  inlet  pressure  gradients,  i.e.  vo  =  poyO-  ~  y),  where  po  =  1, 10,  20.  (a)  The 
initial  distribution  of  QS~  bacteria  shown  in  volume  fractions;  (b)  2D  slice  (z  =  0.5)  of  the 
initial  velocity  profile;  (c)  volume  fractions  of  QS~  bacteria,  QS+  bacteria,  total  bacteria,  as 
well  as  the  concentration  of  AHL,  at  three  selected  inlet  velocities.  The  volume  fraction  of 
QS~  bacteria  and  the  total  bacteria  scales  with  the  inlet  speed,  i.e.,  the  higher  the  inlet  speed 
is,  the  larger  the  volume  fractions  are;  whereas  the  volume  fraction  of  QS+  bacteria  and  the 
AHL  concentration  scales  inversely  with  the  inlet  speed,  i.e.,  the  larger  the  inlet  speed  is,  the 
smaller  the  volume  fraction  and  the  concentration  are. 
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Figure  10:  The  volume  fraction  of  QS+  and  AHL  concentration  at  t  =  15  with  respect  to  three 
selected  inlet  speeds.  This  figure  shows  2D  slices  at  z  =  0.5  for  the  simulation  depicted  in 
Figure  ??.  (a)-(c)  AHL  concentration  distribution  at  t  =  15  for  inlet  velocity  vo  =  Poy(  1  —  y) 
with  po  =  1, 10,20,  respectively;  (d)-(f)  2D  slice  of  the  volume  fraction  of  QS+  (cj)^  )  at 
t  —  15  for  inlet  velocity  vq  =  poy(l  —  y)  with  po  =  1, 10,  20,  respectively. 
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(d)  4>b  i  at  t=25 


(e)  4>b2  at  t=25 


(f)  Nutrient  at  t=25 


Figure  1 1 :  Quorum  sensing  induction  under  hydrodynamic  stress  with  a  weak  nutrient  supply. 
This  figure  shows  quorum  sensing  induction  in  a  flow  cell  with  an  inlet  flow  velocity.  In  this 
simulation,  co  =  0.001,  representing  a  scenario  with  a  weak  nutrient  supply.  Initially  some 
QS~  bacteria  are  attached  to  the  substrate  in  the  middle  of  the  flow  cell,  (a)  The  initial  profile 
of  the  QS~  bacteria;  (b)AHL  concentration  at  x  =  0.5,  y  =  0.0625  and  t=25;  (c)-(f)  2D  slices 
at  z  =  0.5  for  AHL,  ^  1?  ^ 2,  as  well  as  nutrient,  distribution  at  time  t  =  200,  respectively. 
When  nutrient  supplies  are  weak  in  the  entire  flow  field,  QS  induction  is  stronger  downstream 
than  upstream.  Consequently,  there  exist  more  QS+  bacteria  downstream  than  upstream. 
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(a)  c/)bi  at  t=  1.25 


(b)  <frb2  at  t=1.25 


(c)  AHLatt=1.25 


(d)  (pbi  at  t=25 


(e)  4>b2  at  t=25 


(f)  AHL  at  t=25 


(g)  Nutrient  at  t=25 


(h)  EPS  at  t=25 


(i)  v  at  t=25 


Figure  12:  Quorum  sensing  induction  under  hydrodynamic  flow  with  a  strong  nutrient  supply. 
This  figure  shows  quorum  sensing  induction  in  a  flow  cell  with  inlet  velocity  vq  =  10y(l  —  y) 
and  a  strong  nutrient  supply  (co  =  1.0).  In  this  simulation  we  use  the  same  initial  distribution 
of  QS~  as  used  in  Figure  ??  and  characteristic  time  to  =  103  seconds,  (a)-(c)  show  the  2D 
slices  of  the  volume  fraction  of  QS~  bacteria,  QS+  bacteria  and  the  concentration  of  AHL  at 
z  =  0.5  and  t  =  1.25,  respectively,  (d)-(i).  depict  the  volume  fraction  of  QS~  bacteria,  QS+ 
bacteria,  the  concentration  of  AHL,  the  concentration  of  nutrient,  the  volume  fraction  of  EPS, 
and  the  velocity  field  at  z  =  0.5  and  t  =  25,  respectively.  In  this  case,  the  strong  nutrient  supply 
dominates  the  growth  of  bacteria  so  that  QS+  grows  the  fastest  near  the  upstream  location  than 
the  downstream  one. 
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(d)  2D  slice  of  4>n  at  z  =  0.35 


Figure  13:  Hydrodynamic  effects  on  biofilm  morphology  and  the  distribution  of  QS  molecules. 
This  figure  shows  two  grown  up  biofilm  colonies  in  a  shear  flow  in  a  flow  cell,  where  the 
characteristic  time  to  =  1  second,  (a).  The  initial  profile  of  the  two  biofilm  colonies,  (b).  The 
two  colonies  at  t  =  31  in  the  shear  flow.  (c).  Distribution  of  QS  molecules,  AHL.  (d).  A  2D 
slice  of  the  biofilm  colonies  at  z  =  0.35  and  t  —  31.  (e).  A  2D  slice  of  the  AHL  distribution  at 
z  0.35  and  t  31. 
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Figure  14:  An  anatomy  of  the  hydrodynamic  effect  on  biofilm  morphology,  (a)-(c)  show  the 
first  biofilm  colony  in  Figure  (??)  at  time  t  =  0, 18, 42,  respectively.  We  denote  negative  value 
by  red  and  positive  one  by  blue  in  the  rest  of  the  subfigures.  The  data  shown  are  the  solution 
at  t  =  18.  (D)  The  second  normal  stress  difference  Ni  =  (/)n{^yy  —  tzz).  (E)  The  first  normal 
stress  difference  N  +  2  =  (j)n(Txx  —  ryy).  (F)  The  shear  stress  in  xy  plane  (j)nTXy  (G)  The 
shear  stress  in  xz  plane  (f>nTxz •  (H)  The  shear  stress  in  yz  plane  (/)nTyz-  The  first  and  second 
normal  stress  difference  are  of  opposite  signs  in  most  part  of  the  colony.  The  primary  shear 
stress  4>n^xy  is  always  positive  while  the  two  secondary  shear  stress  components  can  take  on 
either  signs. 
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